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D . Abstract 

, In the framework of three generations, we consider the CP violation in neutrino oscillation with 

qq ■ matter effects. At first, we show that the non-perturbative effects of two small parameters, Am^i/Am^ 

and sin #13, become more than 50% in certain ranges of energy and baseline length. This means that 
the non-perturbative effects should be considered in detailed analysis in the long baseline experiments. 
' Next, we propose a method to include these effects in approximate formulas for oscillation probabilities. 

^ Assuming the two natural conditions, 823 = 45° and the fact that the matter density is symmetric, a set 

of approximate formulas, which involve the non-perturbative effects, has been derived in all channels. 

t— I ■ 
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^ ■ 1 Introduction 

. Last year, the direct evidence for neutrino oscillation 1] was found in three kinds of experiments, namely 

atmospheric neutrino experiment [2], reactor neutrino experiment jS] and K2K experiment 4 . In these 
experiments, the dip of neutrino oscillation and the energy dependence of the probability, were observed. 
The possibilities of neutrino decay and neutrino decoherence jS] are excluded by these results and it is 
found that the only solution for the solar and the atmospheric neutrino problem is neutrino oscillation. The 
observation of the dip also means that the neutrino experiments herald in a new era of precise measurements, 
because the effect, which disappears by averaging out the time- varying part on the neutrino energy, has been 
observed in these experiments for the first time [21 01 El- Solar neutrino parameters have been also accurately 
determined by recent neutrino experiments such as SNO and KamLAND |S] . 
' From the results of the past experiments, it was found that the solar neutrino deficit is explained by the 

5^ ! Large Mixing Angle (LMA) MSW solution [71 151 ITU1 [TTj . 

Am^ - 8.0 x KT 5 eV 2 , sin 2 2(9 12 ~ 0.8, (1.1) 

where the mass squared difference is defined by Am^ = rnf — ?n 2 . It was obtained that 

|Am^| - 2.0 x 10" 3 eV 2 , sin 2 29 23 ~ 1.0 (1.2) 

from the atmospheric neutrino experiment |12j . Furthermore, the upper bound of the 1-3 mixing angle, 
sin 6>i3 is given by 

sin 2 26»i3 < 0.2 (1.3) 

from the reactor experiment |13| . The next step for neutrino physics is the determination of sin $13, the 
sign of Am| x and CP phase S. In particular, the measurement of the leptonic CP phase is one of the most 
important themes from the viewpoint of the origin of the universe. CP violation has been investigated also 
in quark sector for the first time and the Kobayashi-Maskawa theory has been established However, it 
has been found that the CP violation in quark sector is too small to generate the sufficient baryon number 
in the universe |15j , because the electroweak symmetry breaking is not the first phase transition as the Higgs 
particle is too heavy. This means that the origin of baryon asymmetry of the universe is not a CP violation 
from the KM phase and an extra source of CP violation is needed. One of the alternatives is the generation of 
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a baryon number due to the leptonic CP violation ■ The possibility of this scenario has been investigated 
by many researchers |17| . 

In order to attain the next step, the long baseline experiments like superbeam experiments |18| and 
neutrino factory experiments are planned. In these experiments, the earth matter effects disturb the 
observation of the CP violation because the matter in the earth is not CP invariant and generate the effects 
of fake CP violation. Therefore, it is very important to understand the earth matter effects in neutrino 
oscillation experiments. 

Here, summarizing the results of the atmospheric, solar and reactor neutrino experiments, there are two 
small parameters 

- 0.04, (1.4) 

< 0.23. (1.5) 

The magnitude of these small parameters is most important for measuring the CP violation, because it 
cannot be observed, if one of these parameters vanishes. As the LMA MSW solution was chosen to explain 
the results of the solar neutrino experiments, a reduced to the largest value compared to other solutions. 
This means that the LMA MSW solution opens the door for measuring the leptonic CP violation. If S13 is 
too small, it will be impossible to observe the CP violation. Therefore the magnitude of S13 controls whether 
the leptonic CP violation can be observed or not. 

Let us briefly review the approximate formulas using the small parameter a or S13 and the related papers. 
At first using the perturbation of oscillation probability in a, the magnitude of the fake CP violation by the 
matter effects has been investigated in [23 EH E2 ESI EH- Furthe rmorc, by expanding the matter potential 
to the Fourier mode, it has been shown in |25l I2fi| that the mode with large wavelength mainly contributes 
to the oscillation probability. Higher order perturbative calculations have been performed by [23 EH1 ■ The 
perturbation in S13 has been investigated in [2H]. The perturbation in both a and S13 has also been studied 
in [301 [3J and this method has been extended to all channels in (32] . 

Next let us review the remarkable features related to the leptonic CP violation. In the case of constant 
matter density, the notable identity JA12A23A31 = JA 12 A 2 3A 31 has been found in [33J [3J] [35] , where J 
is the Jarlskog factor related to the leptonic CP violation, Ay means Am- J -/(2£) and tilde stands for the 
quantities in matter. In addition, it has been pointed out that the oscillation probability in matter almost 
coincide with that in vacuum under the certain condition, which is called vacuum mimicking phenomena, 
and the method to solve the problem on the fake CP violation by using the phenomena is discussed in detail 
|3fil 1371 I3%| . Furthermore, it can be applied to the future long baseline experiments by using the statistical 
method explained by [35] 0U] ■ 

In a previous series of papers ^2 B2I ESI SI C3 we have considered the three generation neutrino 
oscillation in matter and have shown that the CP dependence of the oscillation probabilities are derived 
exactly 0T]. In the case that v e is included in the initial or final state, the CP dependence is given by 



P{v e ^v e ) = C ee , (1.6) 
P{y a — ► up) = A a p cos<5 + B a p sin<5 + C a p, (1.7) 

and in the case that both the initial and final state are v a ,vp = v^,v Tl the CP dependence is given by 

P{v<x — ► vp) = A a p cos 5 + B a p sin S + C a p + D a/ 3 cos 26 + E a p sin 26, (1.8) 



where the coefficients A a p ~ E a p are independent of the CP phase. We have also shown that these coefficients 
can be calculated exactly in constant matter and then the approximate formulas are derived in a simple way 
[421 143j . Furthermore, we proposed a new method for approximating these coefficients in the case of non- 
constant matter density 02], and then applied it to the earth matter |45j . 

In this paper, at first within the framework of two generations, it has been shown that perturbation of the 
small mixing angle is not effective near the MSW resonance point. This means that the non-perturbative 
effects by the small mixing angle is important in the MSW resonance region. Next, we consider non- 
perturbative effects of Am^/Am^ and sin #13 in the three generation neutrino oscillation. The importance 
of the non-perturbative effects is shown by comparing the exact numerical calculation with the perturbative 
expansion of the small parameters. Furthermore, we consider the method for deriving the approximate 
formulas in which the non-perturbative effects are taken into account. In our previous paper |44| . the 
approximate formulas for P{v e — > Vy) have been derived. These formulas are effective for both MSW 



Amii 
Am^ 
S13 = sin 0i3 



2 



resonance regions. However, there is a problem because this method cannot be extended to other channels 
P{ v ji ~ * Vt) an d so on. In order to solve this problem, we assume the two natural conditions, 623 = 45° 
and the symmetric matter potential. Under these conditions, we derive the approximate formulas for all 
channels, including non-perturbative effects of the two small parameters. These formulas are useful to solve 
the problem of parameter degeneracy. - 



2 Non-perturbative Effect by Small Mixing Angle 

In this section, we discuss the perturbative expansion of a small mixing angle in two generation neutrino 
oscillation. Although we discussed the perturbation of small parameters in our previous papers ^J0|)], m 
order to clarify the physical meaning, we consider the perturbation due to a small mixing angle within the 
framework of two generations. Then, we show that the perturbation breaks down in the MSW resonance 
region even if the mixing angle is small. 

2.1 MSW Resonance of Probability in Two Generations 

In this subsection, we consider the two generation neutrino oscillation and we choose the energy region 
and the baseline length in which the MSW resonance occurs. Let us start from the Hamiltonian in constant 
matter 

H = Odiag(0, A)0 T +diag(a,0) (2.1) 
= Odiag(A 1 ,A 2 )0 T , (2.2) 

where the matter potential is defined by a = y/2GpN e . Gf is the Fermi constant and N e is the electron 
density in matter. The matrix O is mixing matrix as 

0=( C0S % sin M, (2.3) 
y — sin 6 cos J v ' 

where A = Am 2 /2E and the quantities with tilde stand for the quantities in matter. Diagonalizing l|2.1|l to 
(|2.2|l . the effective masses Xi(i = 1,2) and effective mixing angle are determined. If we use the notation 
A = A 2 — Ai as the mass squared difference, there is a relation between the mass squared difference and the 
mixing angles as 



''cob 20- sin 2 20. (2.4) 



A sm 20 V V A 
Using these quantities in matter, the oscillation probability is given by 



P = sin 2 20 sin 2 — . (2.5) 

2 v ; 

The oscillating part with L/E of this probability becomes large if the condition 

AL , . 

sin— ~1 (2.6) 

is satisfied. On the other hand, the condition for the maximal effective mixing angle is given by 

sin 20-1. (2.7) 

In the case of small mixing angle, this condition is rewritten as a = A cos 20 — A, and furthermore we define 
the resonance energy by 

Am 2 

E~—. (2.8) 



We also define the resonance length by 



L — . (2.9) 

a sin 
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For the case of sinO = 0.16, which is the upper bound in the CHOOZ experiment, the resonance length is 
roughly estimated as 10000 km. This means that in near future it is impossible to realize the long baseline 
experiments such that the baseline length from beam source to the detector is nearly equal to the resonance 
length. However, it has been shown JSH] that matter effects exist even if the baseline length is shorter than 
the resonance length. Therefore, we use L = 6000 km as the baseline length in the later sections. 

2.2 Perturbation due to Small Mixing Angle 

Next, let us consider the expansion of the effective mass A and the effective mixing angle sin 29 by a 
small mixing angle sin 9. We show that although the effective mass and the effective mixing angle diverge 
in the MSW resonance energy region, the oscillation probability, which is a function of these two quantities, 
converges. 

At first, the effective mass is expanded as 

^ sin 2 g+ o| f A ' sin 4 g + ---. (2.10) 



1 |A-o| 2|A- 

One can see from this result that other terms than the first term diverge. The higher order term have larger 
divergence near the MSW resonance. The effective mixing angle is expanded as 

. ~ A sin 29 ( 2aA . 2 n 3a 2 A 2 . , n \ .„ . 

sin20= — 1 - — -^sin 2 <9 + — s in 4 6> + --- . (2.11) 

A- a V (A -a) 2 2(A - a) 4 / v ; 



The condition 



sin 0<^_^ (2.12) 



2VaA 



is needed for sin 2$ to converge the finite value. However, this condition cannot be satisfied in the MSW 
resonance region defined by A ~ a, even if sin 9 is small. This means that the above perturbation series 
diverges. In the expansion for the effective mass and the effective mixing angle, the coefficients become large, 
even if these quantities are expanded by the small mixing angle. 

Next, let us consider the oscillation probability and let us demonstrate that the oscillation probability 
reaches a finite value, where the divergences due to the effective mass and the effective mixing angle are 
canceled out by each other. Substituting l|2.10[) and (|2.11|) into (|2.5|l . we obtain 



A 2 sin 2 26> 2 (A - a)L 
■ sm — — 



(A - a) 2 
A 2 sin 2 29 
(A - a) 2 



4aAsm 2 6» 2 (A - a)L aALsm 2 9 

■— sin 1 sin(A — a)L 

(A -a) 2 2 A - a 



+ ■■•• (2.13) 



In the limit, A ^ a, it is found that the oscillation probability becomes finite as 



P ~ cos 2 9 ^sin 2 9a 2 L 2 - ^ sin 4 6»a 4 L 4 + •• . (2.14) 

From this equation, the oscillation probability becomes finite and the perturbation is a good approximation 
if 



L<— l —. (2.15) 
asm 9 

As you see from (12.91 . this is the condition that the baseline length is shorter than the resonance length. 

Next, let us investigate the magnitude of non-perturbative effects numerically. We use the following 
parameters, Am 2 = 2.0 x 10 -3 eV 2 and sin9 = 0.16. We set the baseline length, L = 6000 km and the 
energy region, 1 GcV < E < 50 GeV, to include the MSW resonance energy. Furthermore we choose a 
density of p = 4 g/cm 3 . 

At first, in figure 1(a) the level crossing of two eigenvalues is plotted. It is shown that the crossing energy 
is about 6-7 GeV, which corresponds to the MSW resonance energy. Next, in figure 1(b) we compare the 
oscillation probability calculated by perturbation with the one by numerical calculation. These figures show 
that the perturbation breaks down around the MSW resonance energy. The results of this subsection arc 
summarized as 
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(a) Level crossing (b) Probability 




Figure 1: Comparison of the pcrturbative value with the exact one in the two generation neutrino oscillation 
probability. In (a), the energy dependence of two eigenvalues is plotted. In (b), the dotted and solid line 
show the values by the pcrturbative and numerical calculations, respectively. 



1. The perturbativc expansion in the small mixing angle breaks down around the MSW resonance because 
the perturbation because the perturbation series diverges. The coefficients of this expansion become 
larger around the MSW resonance. The divergence included in the effective mass cancels with that in 
the effective mixing angle, and as a result, the value of the oscillation probability reaches a finite value. 
Term of cq. (|2.10|l and (|2.11|) cancel with each other. 

2. Although the divergences of the effective mass and the effective mixing angle in the perturbative ex- 
pansion cancel in the oscillation probability, the finite value of the probability differs from that by 
numerical calculation. The perturbation around the MSW resonance energy becomes a good approx- 
imation, if the baseline length is shorter than the resonance length as seen from l|2.15ll . However, we 
need to take higher order terms of the perturbation into account, when the baseline length is longer, 
namely when non-perturbativc effects become important. 

3 Extension of method to Approximate Oscillation Probabilities 

In this section, we consider the matter effects in three generation neutrino oscillation. At first, we 
review that the 2-3 mixing angle #23 and the CP phase S can be separated from matter effects in the 
oscillation probability [41| . This means that the matter effects appear through the remained four parameters. 
Furthermore, these four parameters can be separated to two set of parameters and each set is related to the 
phenomena in low and high energy as 

(012, Am^) : Low Energy Phenomenon (3-1) 
(#i3, Am^) : High Energy Phenomenon. (3-2) 

This separation means that the parameters for the solar neutrino and those for the atmospheric neutrino 
are almost independent to each other. We propose the method deriving the approximate formulas simply 
by using this feature. 

3.1 Definition of Low and High Energy Regions 

In this subsection, we define the low energy and the high energy Hamiltonians in the small quantity 
limit when S13 or a approximate zero. Although these Hamiltonian have been already introduced in our 
earlier papers [441 145j . we review them here, as they are used in later section. 
It is noted that H (t) satisfies the relation 

H(t) = 23 TH'{t)T^Ol 3 (3.3) 

where H ' is given by 

H' = O 13 O 12 diag(0,A 21 ,A3i)Of 2 Of3 + diag(a(i) ) 0,0). (3.4) 
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This means that the 1-2 and 1-3 mixing angles are separated from the 2-3 mixing and the CP phase, as 
explained in detail in Appendix A. In this Appendix A, we derive the same result as that derived from this 
section from another point of view. Taking the limit S13 — ► 0, the Hamiltonian reduces to the two generation 
Hamiltonian as 

H e = lim H' (3.5) 

Bis-t-Q 

= Oi 2 diag(0,A 21 ,A 31 )Of 2 + diag(a(t),0,0) (3.6) 

A 2 is? 2 + a(0 A21S12C12 \ 
A21S12C12 A 2 ic? 2 (3.7) 
A 3 i / 

This means that the third generation is now separated from the first and the second generation. As seen 
from this Hamiltonian (|3.7|l . the components except for H* T , depend only on (0i 2 ,A2i). We call H l the 
low energy Hamiltonian. On the other hand, taking the limit a — > 0, the Hamiltonian reduces to the two 
generation Hamiltonian as 

H h = lim H' (3.8) 

a— i-0 

= O 13 diag(0,0,A 3 i)O?3+diag( a (t),0,0) (3.9) 

A31S13 + a(t) A 31 si 3 ci 3 \ 

(3.10) 

A31S13C13 A 3 ic^ 3 / 

This means that the second generation is also separated from the two others. This Hamiltonian 13.1l)fl is 
expressed by only the parameters (#13, A31). We call H h high energy Hamiltonian. Next, let us define the 
high and low energy regions described by H h and H e . We first calculate the MSW resonance energy because 
the MSW effect is the most important in matter effects. In the case of L = 6000 km, which we use later, the 
average matter potential is calculated as p = 3.91 g/cm 3 . By using this value, we obtain the high energy 
MSW resonance as E h = Am^/a ~ 5 GeV and the low energy MSW resonance as E = Am^/a ~ 0.2 GcV. 
From these results, we regard E ~ lGeV as the boundary energy of low and high energy regions. Therefore, 
we define the high as E > 1 GeV and the low energy regions as E < 1 GcV. 



3.2 Order Counting of Amplitude on a and s 13 

In this subsection, we investigate how the amplitude S", which is defined by the primed Hamiltonian 
(|3.4|) . depends on the two small parameters a and S13. Before, we have already clarified some general features 
of S' related to the order of a and S13, and the dependences on S13 and a for particular amplitudes S^ e and 
S' Te have been given in our previous papers |441 145| . We investigate now the dependences on S13 and a for 
all amplitudes. 

At first, we represent the explicit form of the Hamiltonian, when the 2-3 mixing angle and the CP phase 
are factored out as 

H'(t) = O 13 Oi2diag(0,A2i,A 3 i)O 1 2O ; r3+diag(a(t),0,0) (3.11) 

(A2ic 13 s? 2 + A3isf 3 + a(t) A21C13S12C12 -A 2 iCi3Si3sf 2 + A31S13C13 \ 
A21C13S12C12 A 2 ic? 2 -A21S13S12C12 . (3.12) 

-A 2 ici 3 si 3 s^ 2 + A31S13C13 -A21S13S12C12 A 2 isl 3 sl 2 + A 31 cj 3 J 

The components of this Hamiltonian depend on a and si 3 as 

/ 0(1) O(a) 0(s 13 ) \ 
H'(t) = O(a) 0{a) 0{as 13 ) . (3.13) 
V O(aia) 0(as 13 ) 0(1) / 

From this result, we can see that non-diagonal components are small compared to the diagonal components. 
We also understand that H'^ is the smallest component and H' H' eT are the next smaller components. We 
should note the salient feature that the result of this order counting holds in H n for arbitrary n. Namely, 
we obtain 

/ 0(1) 0(a) 0(s 13 ) \ 
(H'(t)) n = 0(a) 0(a 2 ) 0(as 13 ) for n=l,2,---. (3.14) 
V O(sia) 0(as 13 ) 0(1) / 
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According to this result, the order of the amplitude S'(t) for two small parameters a and s± 3 is given by 

r f . / 0(1) 0(a) 0(s 13 ) \ 

S"(t) = Tcxp<^ -i / = 0(a) 0(1) 0(as 13 ) . (3.15) 

17 J V 0(s 13 ) 0(as 13 ) 0(1) / 

This result is almost the same as that of the original Hamiltonian. Furthermore, we consider the general 
features derived from the original Hamiltonian. The #13 dependence of this Hamiltonian is described as 

(even even odd \ 
even even odd (3.16) 
odd odd even J 

and this dependence does not change for (H') n , because 

(even even odd \ 
even even odd for n=l,2,---. (3.17) 
odd odd even J 

Due to this result, the amplitude S'(t) has the same structure, 



ft- \ I even even odd 

S' = Texp"^ -i / H'(t)dt > = even even odd 

> \ odd odd even 



(3.18) 



even 

This is a general feature, which holds in arbitrary matter profile. 



3.3 Proposal of Simple Method 

In the previous subsection, we have shown the general features (I3.15f) and (|3.18|) for the amplitude S'(t) 
related to the almost vanishing parameters Si 3 and a. However, we cannot calculate S'(t) by using only this 
features. In this subsection, we propose a generalized method for the calculation. Let us consider if there is 
an approximation available for both region, low and high energy. After expanding the amplitude S' on the 
two small parameters a and s\ 3 , we can arrange this as 

S' = 0(1) + O(a) + 0(si 3 ) + 0(a 2 ) + 0(as 13 ) + 0(s 2 13 ) + ■ ■ ■ (3.19) 

= (O(l) + 0(a) + 0(a 2 ) + •••) + (O(l) + 0(.s 13 ) + 0(8%) + • • • ) 

- 0(l)+0(as 13 )+0(a 2 Sl3 )+0(as 2 13 ) + --- (3.20) 
= limS"+limS"- lim S' + 0(as 13 ) + 0(a 2 s 13 ) + ■ ■ ■ (3.21) 

= S e + S h -S d + 0(as 13 )+0(a 2 s 13 )+0(as 2 13 ) + ---, (3.22) 

where S e , S h and S d arc defined by 

S e = lines' = Texp j -i J H e dtj (3.23) 

(3.24) 



S h = \iinS' = Tcxp^-ijH h dt^ 
S d = lim q S' = diag ^cxp i-i J a(t)dt^ , 1, e 



(3.25) 



respectively. S l (S h ) corresponds to the amplitudes, which gives the main contribution in low (high) energy. 
The term S d counts twice, because contributions to this term comes from both, low energy and high energy 
terms. Therefore, we subtract this contribution and approximate the amplitude as 

S' ~ S* + S h - S d (3.26) 

ignoring higher order terms. Let us discuss this approximation, which is used to derive our main result here. 
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In (|3.19|) - (|3.22|) . the higher order terms in a and S13 are included in S e and S h . The reason for including 
the higher order terms is to take into account non-perturbative effects, which become important in the low 
and high energy MSW resonance region as discussed in section 2. On the other hand, we ignore those higher 
order terms, which are proportional to both a and S13. For example, in the case of second order of the small 
parameters, a and S13, we ignore only the mixed 0(asi3) term among the three terms with second order 
0(a 2 ), 0(s 2 3 ) and 0(asis). This procedure is more appropriate than the usual perturbation, because both 
non-pcrturbativc effects on a small a in the low energy region and on a small S13 in the high energy region 
can be included in our approximation. However, as the derivation of the approximation l|3.26|l is not exact, 
we need to check this later numerically. In the previous subsection, the parity of the matrix elements related 
to S13 has been derived. The equations H3.15fl . I|3.18[l and (|3.26[) lead to the magnitude of the correction for 
the amplitudes as 



^ e = t+°K) (3-27) 

S' Te = S'* e + 0(as 13 ) (3.28) 

S' Tfl = 0(as 13 ) (3.29) 

S' ee = Si e + S' e l e ~Si + 0(asl 3 ) (3.30) 

= S e ^ + 0(as 2 13 ) (3.31) 

S' TT = S? T + 0(as 2 3 ). (3.32) 



If we ignore the higher order terms which arc proportional to both, a and S13, in these equations, we obtain 
approximate formulas by using the two generation amplitudes. The main term for S'^, S'^ is approximated 
by the low-energy amplitude as seen from (|3.27|l and (|3.31|) . On the other hand, the main terms for S' Te and 
S' TT arc approximated by the high-energy amplitude as derived from l|3.28|l and H3.32|) . and these features 
come from eq. (|3.15|) . As seen from ()3.27|) - (|3.32|) . these are expressed by only two generation amplitudes 
and have the advantage of simplicity. The precision of the approximation depends on the values of S13 and 
a. If the value of S13 is smaller than the upper bound derived by the CHOOZ experiment, the precision of 
approximation becomes better. It should be mentioned that the method described in this subsection does 
not need the assumption of constant matter density. 

Next, we show that the results using the approximate formulas <|3.27|) - (I3.32|) are in excellent agreement 
with the numerical calculations. We choose the Preliminary Reference Earth Model (PREM) as an earth 
matter density model and compare the amplitudes in all channels calculated from our approximate formulas 
with the numerical calculation. Here, Am^ = 8.3 x 10~ 5 cV 2 , Am§ x = 2.0 x 10~ 3 eV 2 , sin 2 26*i 2 = 0.8 and 
sin 6*i3 = 0.23 are chosen. Furthermore, we set the baseline length as L = 6000 km, a length, for which 
the MSW effect becomes significant, and the energy region as 1 GcV < E < 20 GeV, for which the MSW 
resonance energy appears. 

We compare our formulas with the numerical calculation in Figure 2. One can see in the following that 
some remarkable features occur. At first, the four amplitudes |5^ e |, | jS"4 e | , \S'„J and \S' TT \ coincide with the 
numerical calculation with a good precision. This happens, because there is no first order correction of 
s 13 from (|3.27|) and H3.30fl - H3.32fl . Next, the low-energy part of \S' Te \ differs from the numerical calculation 
only a little, which can be understood from the cq. (|3.28|) . Furthermore, our approximation for \S' T ^\ is 
not at all in agreement with the numerical calculation. Although the value of this amplitude is exactly 
zero in our approximation as seen from (|3.29() . the actual magnitude of this amplitude attains 0.02 in the 
low energy region from Figure 1. It is noted that this value is almost the same as the value expected 
from the order counting 0(asi 3 ) ~ 0.01. Next, we would like to derive the approximate formulas of the 
oscillation probabilities from the amplitudes obtained here, however, there is a problem. As seen from eqs. 
(IA.32(1 - (|A.49|) in Appendix A, we cannot obtain the approximate formulas for the CP dependence of the 
probabilities Piy^ — > v^^Piy^ — > v r ) and P(v T — ► v T ). The reason is that the CP dependence in these 
channels is directly proportional to S^ T . However there is a method to calculate these indirectly by using 
the unitarity, even if we cannot directly obtain the amplitude S'^ T , as we will show in section 4. 

3.4 Discussion 

In this subsection, let us reconsider the method proposed in the previous subsection in more detail. In 
(I3.19f) - (I3.22|) . we ignored the terms of the order 0{as\ 3 ) for the amplitude S' . The reader probably wonder, 
why we ignore the terms of order 0(as\3) for the amplitude S", but not for other quantities, like for example 
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Figure 2: Comparison of our approximate formulas with the numerical calculation. In these figures the 
absolute value of the amplitudes in all channels is plotted in order to compare our formulas with numerical 
calculation. The solid lines show the approximate probabilities calculated from 13.27fl - 13.32fl and the dashed 
lines show the probabilities in the numerical calculation. 



H' and P. Let us demonstrate the case of using the physical quantity Q. Expanding Q on a and S13, we 
obtain 

Q = Q l + Q h - Q d + 0(as w ) + 0(a 2 s 13 ) + ■■■ (3.33) 

by the same procedure as (|3.19|) - (|3.22|) . If we neglect the higher order terms like 0(asi 3 ), we can approximate 
Q as 

Q~Q e + Q h - Q d . (3.34) 

As in the case of the approximated amplitude defined in the previous subsection, Q e = Q l (9\2 1 A21) is the 
main term in low-energy and Q h = Q h (9i 3 , A31) is the main term in high-energy. Q d is the double counting- 
term. It is a method to be able to take into account non-perturbative effects in both of the two MSW 
resonance regions. In principle, this method is effective whatever we choose for the quantity Q, there is just 
a difference in simplicity and the magnitude of error, as discussed in the following. 
We consider the Hamiltonian H' as Q. Namely, H' can be approximated as 

H' ~ H e + H h - H d , (3.35) 

where the double counting term is given by 

ff d = diag(a(i),0,A 3 i). (3.36) 

There is a problem, because approximation became too simple: The form of the solution for the amplitude 
is given by 



S" - Tcxp j-i j (H* + H h - H d )dtj 



(3.37) 



and we cannot simplify this amplitude without calculation of the commutator of H e and H . Thus, the 
direct application of our method for the Hamiltonian needs other approximations to estimate the amplitude 
and this is not effective from the point of the simplicity. Especially, the amplitudes cannot be calculated 
within the framework of the two generation approximation although the precision of this approximation was 
good. 

Next, let us consider the probability P as the quantity Q. In this case, we can approximate as 

P~P e + P h -P d , (3.38) 
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where P l and P h are given by 



pi{h\ 



Texp j-z J H^dtl 



and P d is the identity matrix. As an example, we consider P(y e 
given by 



(3.39) 



The CP phase 6 dependence is 



A efl cos S + B efl sin 6 + C efl , 



(3.40) 



where the coefficients A e ^ and B efM determine the magnitude of the CP violation. On the other hand, the 
CP violation becomes zero in the limit, a — > or si 3 — > 0, as seen from 



A efl = 0(as 13 ), B efi = 0(as 13 ) 



Namely, we obtain 



A 



A h 



A d 



0, B 



Cfl 



r>h 







(3.41) 



(3.42) 



and therefore we cannot calculate quantities like the CP violation, because it is the effects of three generations 
in this approximation. This result holds for all channels. 

To summarize this subsection, if we choose the Hamiltonian H' as Q, the precision of approximation is 
good, but the calculation is not so simple compared to the exact calculation. If we choose the probability 
P as Q, we cannot calculate three generation effects like CP violation. On the other hand, if we choose the 
amplitude S' as Q, we can calculate the three generation effects like CP violation within the framework of 
two generation approximation. 



4 Approximate Formulas for Oscillation Probabilities 

In this section, we calculate the CP dependent terms from to and so on, not determined by the 
method in the previous section, by using the unitarity. After that, we derive the approximate formulas of 
the oscillation probabilities P(y a — > z/g) in arbitrary matter profile without using S' directly. Namely, we 
derive the approximate formulas in all channels by our new method. 

4.1 Unitarity Relation 

We cannot calculate the amplitude S'^ T in the method introduced in the previous section. The reason 
is that the amplitude S' is a very small quantity, which has an order of 0(as 13 ). As seen from <|A.32|I - 
(JA.49fl in Appendix A, it seems that the approximate formulas, including CP violation, of three channels, 
P(i>n — > Vfj), P(y T — > v T ) and P{v^ — > v T ) cannot be derived without directly calculating the amplitude S'^ T . 
However in this subsection we show, that we can derive these probabilities without directly calculating this 
amplitude, if we assume the two natural conditions, 

S23 - c 23 , (4.1) 
S' aP - S' 0a . (4.2) 

The first condition is supported by the best fit value of atmospheric neutrino experiments |12| and the second 
condition holds in one-dimensional models of the earth matter density like PREM or ak-135f. Accordingly, 
the error due to the difference between these conditions and the real situations is considered to be relatively 
small. We perform the analysis under these two conditions in the following. 
At first, we obtain 

B m = -2Im[(S; M c! 3 + S' TT st 3 y(S' Tfl - 5; T )] C23S23 = (4.3) 
from l|A.34J) and Ij4.2|l in the case of the symmetric matter density. In the same way, we obtain 

B TT = 2lm[(S'^sl 3 + S' TT cl 3 y(S' Tf , - S; T )]c 23S23 = (4.4) 
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from (|A.40|I and (|4.2|) . Furthermore, in the case of the symmetric matter density and the maximal 2-3 mixing 
angle 45°, we also obtain 

A, T = -2Re[(^ - S; T )*(SV4 - ^ T ^ 3 )] C23S2 3 = (4.5) 

from (|A.45|) and (|4.2[) . Let us here consider now, how the oscillation probabilities are derived, which are 
related to the amplitude S' but have not been determined in the previous section,. At first, in the proba- 
bility, 

P(yn ->■ = A w cos 8 + sin 8 + C MM + D MM cos 28 + E m sin 28, (4.6) 

the coefficient proportional to cos 8 can be calculated as 

A m = —A^ - A^ T ~ -A efi ~ -2Re[S^S?J«23CB3 (4.7) 

from l|4.5|l and the unitarity relation. Next, let us turn to the probability P(y T — ► v T ). In the probability, 

P(z/ T — > iy r ) = A rr cos <5 + B TT sin 5 + C rr + D TT cos 2<5 + E TT sin 25, (4-8) 

the coefficient of cos 8 can be calculated as 

A TT = -A Te - A Tfl ~ -A eT ~ 2Re[^;^ e ] S23 c 23 (4.9) 

from l|4.5J) and the unitarity relation. 

Finally, let us calculate the probability P{v^ — > i/ T )- In the probability, 

-P(*V - * ^r) = A^ cos <5 + P Mr sin 8 + C^r + D^ T cos 28 + E^ T sin 28, (4-10) 

the coefficient of sin 8 becomes 

B^r = -B^ e - ~ P e „ ~ 2Im[^:^ e ] S23 c 23 (4.11) 

from l|4.3[) and the unitarity relation. We can derive the probability up to the second order of two small 
parameters by using the unitarity relation although we cannot directly calculate S^ T in the previous method. 
In addition, the coefficients of sin 28 and cos 28, D and E, have an order of 

D = 0(a 2 s 2 13 ), E = 0(a 2 s 2 13 ) (4.12) 

in these three channels as derived from l|A.36|) . I|A.37(I . (|A.42|I . (|A.43|) . I|A.48|) and JA.49I and are expected to 
be small. Actually, these coefficients have the second order of 5^ T , and the values are about (0.02) 2 ~ 0.0004 
from Figure 1 in the high energy region related with long baseline experiments. So we ignore them in the 
following section. 



4.2 Approximate Formulas in All Channels 

In this subsection, we present the approximate formulas which are useful in arbitrary matter density 
profile. Ignoring the higher order terms of a and si 3 than the second order, we can present the oscillation 
probabilities for all channels with the amplitudes calculated in two generations. 

At first, let us derive the approximate formulas for P{v e — > v^) and P(v e — > v T ). The approximate 
formula for P(v e — ► v^) has already been derived in our previous paper |44| . We only have to replace the 
amplitudes £' and S' Te in three generations into S^ e and S^ e in two generations. From (|A.24|) - (|A.31J| and 

we obtain 

P(^ e — > i^) = A e)1 cos 8 + B e ^ sin 8 + C ejl (4-13) 

Aep, * 2Re[S%S , ; e \c 23 s 23 , (4.14) 

B efi ~ -2lm[S%S'; e }c 23 s 23 , (4.15) 

Ce, * \SU 2 c 2 23 + \S h Te \ 2 s 2 23 , (4.16) 

P{y e — > v T ) = A eT cos 8 + B eT sin 8 + C er (4-17) 

A e r -2Re[^* e ^ e ]c 23 s 23j (4.18) 

B eT ~ 2Im[^* e ^ e ]c 23S23 , (4.19) 

Cer * \SU 2 4 3 + \S^fc 2 23 ,. (4.20) 
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Eqs. H4.13f) - I|4.16f) are the same as those derived in our previous paper 23 ■ Next, let us derive the approxi- 
mate formulas for P(v e — > v e ). Using (|3.30[) directly, we obtain 

P(v e ^v e ) = C ee = \S'J 2 (4.21) 

* \Sl + S h ee -S d e f. (4.22) 



On the other hand, we obtain 



P{v e ^v e ) = C ee = l-C ell -C eT (4.23) 



eft 

Kf-\S h Te \\ (4.24) 



by using the unitarity relation. This is a different approximate formula than (|4.22() . Thus, there are two 
kinds of expressions (|4.22(l and 1)4. 24|) for P(v e — > v e ). We checked numerically that the expression (|4.24(l 
has a better precision than the expression (|4.22(l . Furthermore, the expression l|4.24(l easily reproduces the 
approximate formula derived with double expansion up to the second order of two small parameters in rcf. 
[221 (second order formula). So we use the expression (|4.24|) in the following. 

Next, let us derive the approximate formula for Piy^ — > v r ). At first we calculate the terms independent 
of the CP phase in this calculation. We can approximate 

Cpr = \S'» T \ 2 sl 3 + 1^/4 + \S'^ - S' TT \ 2 c 2 3 s 2 3 ~ \S^ - S h TT \ 2 c 2 23 s 2 23 (4.25) 

from <|A~47|) and IpOTjl - lpO^ . where we ignore the terms proportional to |<S^ T | 2 = 0(a 2 s\ 3 ). This leads to 
the approximated probability as 

P(v m ->Vt) = B^sinS + C^ (4.26) 

Bfir * 2Im[^;^ e ]c 23S2 3, (4.27) 

C, T ~ \S^ - S? T \ 2 4 S s 2 23 . (4.28) 

Next, let us derive the approximate formulas for P{v^ — > v^) and P{v T — > v T ). From l|A.35|l and H3.31fl - l|3.32[l . 
we obtain 

= l^4+^ T s 2 23 | 2 + (l^| 2 + l^| 2 )4s 2 23 (4-29) 

— l^/i c 23 



1 S h TT s 2 ,\\ (4.30) 



where we neglect the terms proportional to |S^ T | 2 = 0(a 2 s\ 3 ). On the other hand, we obtain another 
expression by using the unitarity relation as 



Cfj,^ — 1 C^ e Cfj, T (4.31) 

Vel u 23 ~ I'-Vel °23 — 



1 - \SL\ 2 ^ \S h T f^ - \S' - S h TT \ 2 c 2 3 s 2 3 (4.32) 



This seems to be different from (14.3011 at a glance, but we confirmed that l|4.30[) and (|4.32|) are the same 
expression by using the unitarity relation. In the following, we use the expression ()4.32|) for the reason that 
this easily reproduces the second order formula and we can check the unitarity. In the same way, C rr is 
given by 

C TT = 1 — C eT — ~ 1 — |S^ e | 2 S23 — |*S> e | 2 C2 3 — — S TT \ 2 C^ 3 s\ 3 (4.33) 

from the unitarity relation. From the result obtained in subsection 4.1, the approximate formulas for P(y^ — > 
v^) and P(y T — > v T ) arc given by 

P{vp -> v») = cos 5 + C MM (4.34) 
~ -2Re[S%S* e }c 23 s 23 , (4.35) 

— ^ — l^e| 2c 23 ~ I^Tel 2 ^ — l^/i/i — ^tt\ 2 c 23 s 23 (4.36) 

P(y T -> v T ) = A^cosJ + CW (4.37) 
A TT ~ 2Re[^* e ^ e ]c 2 3S23, (4.38) 

Crr * l-|^ e | 2 S | 3 -KM 2 4-|^-^ T | 2 ci 3 43- (4-39) 
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These results are one of the main results of this paper. In all channels, we can present the probabilities 
including the CP violation by using the amplitudes calculated in two generations. It is noted that the 
CP violating terms due to the existence of three generations can be calculated from the two generation 
amplitudes. 

Next, let us compare the approximate formulas H4.13M4.20jl . Q4.24fl . |4.2fiM 4 -^ and lj4.34M4.39jl with 
the numerical calculations. We take the PREM (Preliminary Reference Earth Model) as the earth matter 
density profile and compare the approximated values of all probabilities with those calculated numerically. 
We use the same parameters as those used in fig. 1 and sin 2# 2 3 = 1, S = 90°. We set the baseline length, 
L = 6000 km and the energy region, 1 GeV < E < 20 GeV, to include the high energy MSW resonance. 



P{v e -> Vp) 



P(Vfj, — > v T ) 
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1 h U , 


0.8 








. 6 








0.4 








0.2 









1.5 2 3 5 7 10 15 20 





1.5 2 3 5 7 10 15 20 



1.5 2 3 5 7 10 15 20 



1.5 2 3 5 7 10 15 20 



Figure 3: Comparison of our approximate formulas with the numerical calculation. In these figures P(v a — * 
i/g) is plotted in order to compare our approximate formulas with the numerical calculation. The solid lines 
show the approximate probabilities and the dashed lines show the numerical calculation of probabilities. 



We compare our approximate formulas with the numerical calculation in figure 3. One can see some 
features from this figure. The approximated value of probabilities for P(y e — > u^), P(y e — > v T ) and P(y e — ► 
v e ) coincide to the numerical values very well, on the other hand, the remaining three channels of probabilities 
P( u fj, ~ * v t)i P{ v h ~ * v n) an d P(vt Vt) show a small difference between the approximate and the numerical 
value. However, the difference is not so large as in figure 2 and as a first step the result is sufficiently accurate. 



5 Comparison of Our Results with Second Order Formulas 

In this section, we concretely calculate the amplitudes by using the approximate formulas derived in the 
previous section for the case of constant matter and show that simple approximate formulas can be obtained. 
Finally, we demonstrate that the approximate formula derived with double expansion up to the second order 
of the two small parameters (second order formulas) are largely different from the exact values in the MSW 
resonance region under the condition that the baseline length is longer than the oscillation length. 



5.1 Approximate Formulas for Amplitudes 

In the previous section, we have given a method for approximation of the probabilities in three generations 
by amplitudes in two generations. In this subsection, we use the constant matter density profile in order to 
compare our method with other method and investigate the non-perturbative effects. As seen from (j4.13jl - 
(j4.2()jl . Ij4.24jl . lj4.26M 4 -28jl and (j4.34M4.39jl . we only have to calculate four kinds of amplitudes, namely 
^pei Syy, S Te and S TT . 

The low-energy approximate formulas are obtained by taking the limit S13 = and from 

H l = O 12 diag(0,A 21 ,A 31 )Of 2 +diag(a,0,0) (5.1) 
= 0{ 2 di ag (\l\i,A 31 )(0{ 2 ) T . (5.2) 
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The effective masses = 1, 2) and the effective mixing angle 9{ 2 are determined by the diagonalization of 
(|5.1|l to Ij5.2fl . If we define the mass squared difference in matter as A 21 = X 2 — X\ , we obtain the relation 



A e 21 sin26»i2 f no a 



cos20 12 -— + S i^20 12 . (5.3) 
A 2 i sin20j 2 y V A 2 i/ 

Therefore, the amplitude is calculated as 

SL = -i s in29{ 2S m^e X p(-i^±^L] (5.4) 



Si, = (cos^-icos2^ 2S in^)exp(-i^L±^i) (5.5) 

by substituting l|5.2|l into Ij3.23j) . On the other hand, the approximate formulas in high energy are obtained 
by taking the limit a = and we get 

H h = O 13 diag(0,0,A 31 )Of 3 +diag( a ,0,0) (5.6) 
= O^ 3 diag(A' 1 l ,0,A 3 l )(O^) T . (5.7) 

The effective masses (i = 1,3) and the effective mixing angle 6^3 are determined by the diagonalization of 
(|5.6J) to (|5.7J) . If we define the mass squared difference in matter as A 31 = A 3 — A'/, we obtain the relation 



^ = ^ = \/ (cos 20 13 - ^-)\s^29 13 . (5.8) 
A31 sin26>y 3 y V A31/ 

Accordingly, the amplitude can be calculated by substituting (|5.7(l into (j3.24j) as 

S* - -isin2^ 3 sin^|^ex P ('-i^L±^^ (5.9) 

S * t = f C0S ^-i C082 e^sm^y^>(-i^^Ly (5.10) 

As seen from Ij5.4|) and l|5.9|) . and Sy e have simple forms, but the expressions of and S% T are more 
complex than Ij5.5|l and l|5.10|l . 

5.2 Approximate Formulas for Probabilities 

In this subsection, we derive the approximate formulas of the oscillation probabilities in constant matter 
by using the result of the previous section. 

At first, let us consider the case of including electron neutrino in the initial or final state. In this case, 
the probability for any channel can be calculated almost in the same way. The probability P(y e — > v e ) is 
obtained by substituting Ij5.4|) and Ij5.9jl into (|4.24|) as 

P{u e -^v e ) = l- sin 2 26[ 2 sin 2 ^ - sin 2 20% sin 2 (5.11) 

The probability P(v e — > v^) is obtained by substituting (|5.4jl and (|5.9jl into H4.14jH4.16j) as 

P^e—tv^) = A efl cos 5 + B efl sin S + C eM (5-12) 

A^ L A' 1 L A L 
A ell ~ sin 2^ 2 sin 26» 23 sin 2^3 sin sin— |— cos — |— (5.13) 

B^ ~ sin 2^ 2 sin 26*23 sin 2^ 3 sin ^2A^ sin sin (5.14) 

A^ i A ft i 

C efi ~ c 23 sin 2 26*^2 sin 2 — |— + sin 2 26^ 3 sin 2 (5.15) 

The remaining probability P{v e — > v T ) can be calculated from the unitarity relation. Next, let us calculate 
the probabilities for the case, that not all electron neutrinos in the initial and final state are included. Also 
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in this case, the probability for any channel can be calculated almost in the same way. Accordingly, as an 
example, we calculate the probability for muon neutrino to tau neutrino, 



B„ T sin 5 + C^r 



BfiT 

Cut 



ALL . AS, L 



A, 2 L 



sin 26{ 2 sin 26> 23 sin 26^ 3 sin ^1±L s i n t^k^ 



q( _ oh |2 2 : 



23 ■ 



(5.16) 
(5.17) 
(5.18) 



At first, we use the relations, cos2#f 2 = 2cos 2 #f 2 — 1 an d cos26^ 3 = 1 — 2 sin 2 0^ 3 , and we rewrite S^L and 

S clS 



s 



I'll 



nil 



exp i 



exp 



,^21 



L\—2i cos 



2 a£ 
7 12 



sin 



A* 21 L 



exp — z 



.A 21 



; A>\Z 



A' 1 L 

2i sin 2 e\\ sin 



exp — z 



2 

,A 31 



-L 



Then, arranging C MT in the order of the effective mixing angles cos 9\ 2 and sin (9^ 3 , we obtain 

(Ai + A§! + A 32 )L 



r 1 



= sin 2# 93 sin 



C 



2u 



(A £ 21 + A§ 1 + A 32 )i (A| x - A§! - A 32 )Z . A^L 



C 



2f» 



-2 sin 2 26» 23 cos 2 6»f 2 sin 

?sm 2 9tf , in 2 M , in ( A 21 + + A 32 )£ (Aj - A 3 \ 
-z sm z(79:i sm (7i o sm cos 



A 32 )L . A 31 L 
sm — ^ — 



CL = sin 2 26*23 cos 4 l 12 sin 2 



A £ 21 i 



/IT 



2 sin 2# 23 cos 



2 nl 



'12 



sin 



2 M 
7 13 



sm 



sin 2 26» 23 sin 4 6^ sin 2 31 



A^L . A 31 L A 32 L 



sm 



cos ■ 



(5.19) 
(5.20) 

(5.21) 
(5.22) 
(5.23) 

(5.24) 



As we show in the following section, the reason of arranging the terms like this is, because the second order 
formulas can be easily derived. In order to derive the second order formulas, it is sufficient to use C* T , C 2 " 
and C 2 ^.. We can also calculate the other channels Piy^ — > v^) and P{y T — * v T ) in the same way. In a recent 
study, it was found that the channels P(y^ — > v^) and P(y T — > v T ) are largely affected by the earth matter 
in the long baseline [4T1 El9"| . 

We can see from these expressions that the approximate formulas are rather complex for the case not 
including electron neutrino in the initial or final state. Wc also understand from these formulas how matter 
affects the probabilities. Thus, the formulas are expected to be useful for studying matter effects. 



5.3 Large Non-perturbative Effects of small a and s 
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In this subsection, we compare the approximate formulas obtained in the previous subsection with the 
second order formulas numerically and it is shown that the latter have a large difference from the numerical 
value in the MSW resonance region. 

The second order formulas are approximated by the main terms of the expansion and are widely used by 
many authors for their simplicity. In refs. |3L>I 131 j . the formula for P(y e — > i/„) has been derived and later 
on all probabilities were presented in ref. |32j . As examples, the probabilities, P(y e — > v^) and P (ia, — > v T ), 
are taken. For the other channels of probabilities, similar expressions have been obtained. In all channels 
similar results were obtained from comparison with numerical calculations. The second order formula for 
P{y e — > Vfj) is given by 



P(y e 



A, 



A e)1 cos 5 + B efl sin S + C en , 



e 4 l 



B 



as 13 sin 29i 2 sin 26*23 



asi3 sin 26i 2 sin 26> 23 



9A 2 



aL . (A 3i -a)i A 32 L 

, sm — sm cos 

a(A 3 i-a) 2 2 2 



2A 2 ! . aL . (A 31 - 
—n r sm — - sm — 

a(A 31 - a) 2 2 



a)L . A 32 L 
sm — - — 



a 2 c 23 sin 2 26»i 2 ^- sin 2 ^ 
a z 2 



2 2 
S 13 S 23 



4A 2 , 



(A 31 -a)Z 



(A 31 - af 



(5.25) 
(5.26) 

(5.27) 

(5.28) 
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Comparing our approximate formulas l|5.13fl - l|5.15fl with the second order formulas (f5.26f) - (f5.28f) . each term 
corresponds one by one. Actually, the second order formulas <|5.26f) - H5.28fl are derived by expanding our 
approximate formulas (|5.13() - (|5.15() up to the second order in a and s\ 3 |44| . Next, the second order formula 
for P(y^ — > v r ) which has been already derived in ref. |32| is 



An T 



= A^ T cos 5 + B^ T sin 8 + C^ T 
~ asi3 sin 2 26*23 sin 26*12 cos 26*23 



2A 



31 



. A31L 



A 31 - a 



a . A 31 L A 3 i . aL (A 31 - a)L 

— — sm sin — cos 

A31 2 a 2 2 



asi3 sin 26*12 sin 26*23 



2A 2 31 aL . (A 31 -a)L . A 32 i 



, sm — sm ■ 
a(A 31 - a) 2 



■ sm ■ 



(5.29) 

(5.30) 
(5.31) 



and C^t is given by 



■ 2 na -2 A31L 

sm 26*23 sm 



a sin 2 26*23 cos 2 6*12 ( — — — J sin A31L + a 2 sin 2 26*23 cos 4 6*12 ( — — — I cos A31L 



A31 



a z sin 2 26*23 sin 2 26* 12 ( 



A 31 L (A 31 -a)L . aL f A 31 \ A 31 L . . 
sm __ cos _ Sln _ |_ j _ __ sm(A 31 L) 



s 2 3 sin 2 26* 



2:i 



2A 3 1 

A31 - a 



A 31 L aL . (A31 -a)L( A 31 \ aL 
sm __ cos _ sm _ f __ j _ _ sm{ A 31 L) 



(5.32) 



In the next section we show that this formula l|5.32|l can be also derived from our formulas H5.21fl - H5.24fl . 
It is noted that the second order formula (|5.32|) for C^ T is rather complex. Furthermore comparing our 
approximate formula l|5.17|) - (|5.24|) with the second order formula H5.30fH5.32f> . we see that A^ T is not 
included in our formula. The reason is, that A^ T = in the case of maximal mixing angle sin 26*23 = 1 and 
there is no way of calculating this by the method described in this paper. If we consider cos 26*23 as a small 
parameter like a and si 3 , this A^ T has the magnitude of 0{as\ 3 cos 26*23). Therefore, A^ T is proportional 
to the third order of small parameters and is expected to be neglectable. This means that our formula is 
not largely affected by the error due to this term, which cannot be derived from our method. However, as 
this error affects the precision measurement of sin 6*23 by the atmospheric neutrino experiments in future, 
the improvement of this point is important future work. The formulas for the other channels are given in 
ref. |32| . The second order formulas are effective under the following two conditions. 
The first one is for the neutrino energy and is given by 



E > 0.45 GeV 



A? 



'21 



10- 4 eV 2 J 



3 g/cnr 3 



The second one is for the baseline length and is given by 

L« 8000 km (JL) 



10" 4 eV^ 
Am| x 



(5.33) 



(5.34) 



These conditions come from the utilization of perturbative expansion on the two small parameters. The 
detailed discussion arc given in |31j . These approximate formulas are used for the purpose of understanding 
of the results obtained by numerical calculations |501 151] . However, as shown in the next figure, these 
formulas have large difference from the true value in the MSW resonance region, which is considered to be 
the most important region. 

Next, let us compare our formulas H5.12fi - I|5.23fl with the second order formulas (|5.25fl - (f5.32f) in all chan- 
nels by numerical calculation. In order to see the magnitude of the error, we also compare two kinds of 
approximate formulas with the exact values. We set the baseline length as L = 6000 km, where the MSW 
effect becomes large, and the energy region as 1 GeV < E < 20 GeV, where the MSW resonance energy 
is included. Furthermore, the second order formulas are derived only in the case of constant matter, so we 
choose the average density p = 3.91 g/cm 3 of the earth calculated by the PREM. Note that two conditions 
(|5.33fl and (|5.34fi are satisfied in these region. 



16 




1 1.5 2 3 5 7 10 15 20 1 1.5 2 3 5 7 10 15 20 1 1.5 2 3 5 7 10 15 20 

E GeV E GeV E GeV 



P(y e — > t/ e ) P(i/ p -> P(^ r -> i/ T ) 




1 1.5 2 3 5 7 10 15 20 1 1.5 2 3 5 7 10 15 20 1 1.5 2 3 5 7 10 15 20 

E GeV E GeV E GeV 



Figure 4: Comparison of our approximate formulas with the second order formulas and in addition with 
the numerical calculations. The solid, dashed and dotted lines show the probabilities in our approximate 
formulas, those in the numerical calculation and those in the second order formulas, respectively. 

We compare the probabilities calculated by our approximate formulas in all channels with those by the 
second order formulas in addition to numerical calculation in figure 4. One can see the following points from 
this figure. The second order formulas show large differences from the numerical values around 5 GeV, where 
the high energy MSW resonance occurs. In other energy regions they are in good coincidence. The value of 
P(v e — ► v e ) has the largest difference, the probabilities P(v e — > f») and P(y e — ► v T ) have the next largest 
difference, the values of P(i / (U — * f M ) and P(y T — ► v T ) have also significantly large difference, but only the 
probability P(i / (U — * f T ) has a small difference. In addition, these figures show that the difference between 
the second order formulas and the numerical calculation exists even in the two applicable regions 1(5.33(1 and 
((5.34() . Although we do not show a figure, the difference between our approximate formulas (|5.1 2(1 - (|5. 23(1 
and the second order formulas 1(5.25(1 - 1(5.32(1 become more clear out of the two applicable regions 1(5.33(1 and 
((5.34|) . The reason is that our approximate formulas ((5.12|) - ((5.23(l arc applicable even for the case that the 
condition ((5.33|) or ((5.34|> for energy and baseline length does not hold, as confirmed from the comparison 
with the exact numerical calculation. However, the second order formulas are good approximations, when 
the neutrino energy is not near the resonance energy, even if the baseline length is long. 

6 Non-perturbative Effects of Small Parameters Am^/Am^ and 
sin 6 iz 

In this section, we investigate the reason for the difference between the second order formula, which 
contains the approximation with double expansion up to the second order of two small parameters, and the 
numerical calculation around the MSW resonance region as explained in the previous subsection. 
We discuss the non-perturbative effects of small mixing angle more detailed than in section 2. 

6.1 Derivation of the Second Order Formulas 

In this subsection, we investigate how the second order formulas are approximated expanding on a = 
Am^/Arrijj In the previous paper, we have discussed the probability P(v e — > f„), so we calculate the second 
order formula for P(y^ — > u T ) here. The method of calculation is basically the same but the calculation 
itself becomes slightly complex, because we need to calculate the effective mass and the effective mixing 
angle up to the second order of a and S13 in the case of Piv^ — > v T ). In this point, the calculation is not 
straightforward compared with that of P(v e — ► v^) but the method of approximation is the same. Note 
that CL. in 1(5.21(1 does not include the effective mixing angle. For this reason, it is sufficient to expand the 
effective mixing angle up to the zeroth order, but we need to expand the effective mass up to the second 
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order to calculate the probability up to the second order in a and S13, The effective mixing angle is expanded 
up to the zeroth order as 

cos6»i 2 ~ -sin20f, ~ asin20i 2 — (6.1) 
2 2a 



nd h 13 ~ i s i n 2^3~ Sl3 -^- (6.2) 
Z A31 — a 

and the effective mass is expanded up to the second order as 



A 2 

A|i ^ a-acos26»i 2 A 3 i +a 2 sin 2 26»i2^ i (6.3) 

ACL 



A 3 \ ~ A 31 -a + S 2 3 -^i^. (6.4) 
A31 - a 



Here, we should emphasize the following points. Eqs. (|6.1|l and l|6.3[) obtained by the expansion in a 
diverge in the vacuum limit a — ► and eqs. (|6.2|l and l|6.4|) obtained by the expansion in S13 diverge in the 
high energy MSW resonance limit a — > A31. As shown in the following, these divergences cancel and the 
probability has a finite value. Expanding C* in l|5.21J) up to the second order, we obtain 

= sin 2 20 23 sin 2 ^ (6.6) 
= -a sin 2 20 2 3 cos 2 12 f^^\ s i n A 3 iL (6.7) 

= a 2 sin 2 2023 cos 4 12 f^f^) cosA 31 L (6.8) 

C# - a 2 sin 2 2023 sin 2 20 12 sin A 3i L (6.9) 

= s 2 13 sin 2 2023 f - |) sin A31L. (6.10) 

We also expand C 2 ° in (|5^ and C 2 ^ in tfZZfy as 

* -2. 2 sin 2 2023 sin 2 20 12 (^) ' sin ^ cos (Asl ~ sin ^ (6.11) 

C* * ~2 S 2 3 sin 2 2023 (5^) 2 sin ^ cos ^ sin ~ ° )L . (6.12) 

Finally, we obtain (|5.32|l arranging these result order by order. 

Here, let us consider the applicable region of the second order formulas. C^f diverges in the limit a — ► 
and Cj^ diverges in the limit a — ► A31. C 2 ° also diverges in the limit a — > and C 2 * diverges in the limit 
a — > A31. The divergences in a — > and in a — > A31 come from the expansion of the effective masses <|6.3|) 
and (|6.4(l respectively. It seems that the second order formulas do not reduce to those in vacuum due to the 
divergence in a — > and furthermore do not reduce to those in the high energy MSW resonance point due 
to the divergence in a — * A31. However, when we consider the pair 



Cli + Cl a T = -a 2 sin 2 2023 sin 2 20 12 ^1 
^ ^ la 



A 31 L (Asi-o)L . ah /A 31 \ A 31 L . 
sin — ^— cos sm—\— — sm(A 31 L) 



the divergence in a — > cancel and the value converges. The obtained finite value is given by 



(6.13) 



Hm(C^ + CIV) = -o? sin 2 20 23 sin 2 20 12 i sin 2 (6.14) 

Before we expand, and C 2 T have finite values in the limit a — *■ and a — > A31. However, the divergence 
appears in expansion of a and S13. The cancellation of these divergences occurs between Cj^ and C 2 °. This 
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means that the cancellation occurs between the different terms and result in finite values, respectively, at 
first, which is an interesting result. Considering the pair as 



Cl 



C 2b 

(AT ' 



-s? 3 sin 2 26*23 



2A 3 i 
A 3 i -1 



A 31 L aL . (Asi-o)L 

sm cos — sin 

2 2 2 



A31 



A31 



aL 

— sm(A 31 L) 



the divergence in the limit a — ► A31 cancels and the value converges. The finite value is given by 



lim (C£ + C£) = ^3 sin' 2023 



A31L 



(A31L) sin : 



A31L 



sin(A 3 iL) 



(6.15) 



(6.16) 



The cancellation of these divergences occurs between the different terms C^ T and C 2 *, which is also a 
remarkable result. 

We have shown that the second order formulas have finite values in the limit a — > and a — > A31, 
but it is not always the same as that in the numerical calculation. Actually, the difference in fig. 4 in the 
limit a — > A31 , shows that the second order formulas have finite values but they are not in accordance with 
those in the numerical calculation. In order to study this, we compare the three quantities, the numerical 
calculation, our formulas and the second order formulas. We can learn the differences mainly in the vacuum 
limit a — > and the high energy MSW resonance limit a — ► A31 from the comparison. 

At first, let us consider the vacuum limit a — > 0. Furthermore, to simplify the discussion, we consider 
the case of S13 — > 0. The second order formulas in the limits a — > and S13 — > are given by 



lim C^ble) 



2 2 ^^31-^-' 9 9 

sin 26*23 shi a sin 26*23 cos 1 



'12 



sin A31L 



a 2 sin 2 26*23 cos 4 6*i 



A31L 



COS A31L 



a 2 sin 2 26*23 sin 2 26»i 2; 



1 /A31L 



sm 



A31L 



(6.17) 



Next, taking the limit a — ► and S13 — > in our formulas, we obtain 



lim C (exact) 

a,si3^U 



sin 26*23 sin" 



A31L 



2 sin 2 26*23 cos 2 6*12 sin 



A31L (A 2 i - A 3 i)L . A21L 



(6.18) 



2 2 2 

Expanding the oscillating part of H6.18(l in our formula, it leads to (|f> . 1 T|> obtained from the second order 
formula. The condition for the expansion on the oscillating part for sufficiently good approximation is 

2 



L < 



A 2 i ' 



(6.19) 



Next, let us consider the high energy MSW resonance limit a — ► A 31 . In order to simplify the discussion, 
we take the high energy MSW resonance limit a — > A 31 under the condition a — > 0. In the high energy 
MSW resonance limit of the second order formula, we obtain 



lim (-(double) 



a — >A3i ,a — >0 



. 2<ifl -2 A31L 

sin 2023 sm — - — 

2 • 2 oa ( ^3iL 
s 13 sm 26*23 



(A31L) sin : 



A 31 i 



sin(A 3 iL) 



(6.20) 



Next, taking the limit a — ► A 31 and a — > in our formulas, we obtain 



lim C, ( , 



a — ^A 3 i ,a — >0 



■ 2 oa ■ 2 (1 + Sl3)A 3 iL 

sm 26* 23 sm 



sin 2 26* 23 (1 - s 2 3 ) sin 



(l + si 3 )A 31 L (l-si 3 )A 3 iL 



■ cos 



sin(s 13 A 31 L). (6.21) 



2 4 

By expanding the oscillating part obtained from our formula l|6.21|l . it is shown that this coincides with that 
from the second order formula l|6.20|l . The condition for the expansion of the oscillating part for a sufficient 
approximation is given by 

2 2 

L < — = (6.22) 

si 3 A 3i si 3 a 
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If the baseline length is shorter than that obtained from above condition, the second order formula becomes 
a good approximation. We obtain the following results about the perturbative expansion on the small 
parameters a and S13. 

1. The perturbative expansion in a actually corresponds to the expansion in A21/CI. This constrains the 
applicable energy for the approximate formulas. If we expand in the parameter A21/C1, the effective 
mass Afi and the effective mixing angle sin2(9^ 2 diverge in the vacuum limit a —> 0. However, these 
divergences cancel out each other in the calculation of the oscillation probability. Thus, the probability 
has a finite value, but the value largely differs from the numerical calculation in low-energy. The 
magnitude of this difference becomes large and serious in the case of small mixing angles and in 
low-energy long baseline experiments. 

2. If we expand in the small mixing angle S13, the effective mass A^ and the effective mixing angle 
sin 28 1 3 diverge in the MSW resonance energy limit a — > A31. However, these divergences also cancel 
each other out in the calculation of the oscillation probability. Thus, the probability has a finite value, 
but the value largely differs from the numerical calculation in the high-energy MSW resonance region. 
This means that the second order formulas cannot be used in the high energy MSW resonance region. 

In two generations, we can calculate the oscillation probabilities exactly by solving the second order equation. 
So, we do not need the perturbative expansion. On the other hand, the construction of the approximate 
formulas applicable to arbitrary matter density profile is very difficult in three generations. Therefore, we 
need to expand on the small parameters a and S13. 

6.2 Discussion 

We have shown that the double expansion formulas up to the second order in the two small parameters a 
and S13 does not give a good approximation in the MSW resonance region. This is because the coefficients of 
the small parameters have large values in the MSW resonance region. In this subsection, let us discuss some 
methods proposed up to present to solve this problem. The Hamiltonian H' is written by four parameters. 
The two parameters (Am^,^) control the physics mainly in the low-energy region and the other two 
parameters (Am^, 9 13 ) control the physics mainly in the high-energy region. In other words, the magnitude 
of a determines low-energy phenomena and the magnitude of S13 determines high-energy phenomena. Both 
of these parameters are very small but the energy region, where the expansion converges, is different. This 
means that we need to treat the applicable energy region carefully when we expand on these two parameters. 
There are several methods in order to take into account the higher order terms of a and S13 for example 

1. exact formulas in constant matter density profile 

2. reduction formulas taking into account the two generation part exactly 

In the first method, there does not exist any error generated from the perturbative expansion, because of 
the exact treatment of both a and S13 |42j . Furthermore, non-pcrturbative effects can be easily investigated 
by using these exact formulas. The second method was introduced in our previous paper |44j . In this 
method, we try to include the higher order terms of a and S13 partially, except for the terms including the 
product of two small parameters. This method includes the higher order terms of a and S13 and is simply 
and applicable even in the case of arbitrary matter density profile [441 14~5] . Although this method uses only 
the second order approximation of the amplitude, it has the notable feature that the third order (three 
generation) effects such as CP violation can be calculated. 

7 Summary 

In this paper, we consider the method how to approximates the neutrino oscillation probabilities in 
matter under three generations and the obtained results are summarized as following. 

1. In the framework of two generation neutrino oscillation, we discuss the applicable region of the per- 
turbative expansion on the small mixing angle in matter. The result of the perturbation differs largely 
from the exact numerical calculation in the MSW resonance point. This means that non-perturbative 
effects arc important even for the neutrino oscillation in two generations. 
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2. We extend the method to calculate the approximate formulas, in which non-perturbative effects 
of the small parameters Am^/Amjj and sin #13, to all channels. Under the conditions, #23 = 45° and 
the symmetric matter density profile, we derive simple approximate formulas of the probabilities in all 
channels by using the unitary relation. Although all these approximate formulas are expressed by the 
amplitudes calculated within the framework of two generations, it has a notable feature that the three 
generation effects such as CP violation can also be calculated. 

3. In the three generation neutrino oscillation with matter, we investigate non-perturbative effects of the 
two small parameters Am^/Am 2 ^ and sin #13. We compare our approximate formulas with those from 
the double expansion, which include the terms up to the second order in the low and high energy MSW 
resonance regions. The obtained result is that the second order formulas show large differences from 
the exact numerical calculation, which means that non-perturbative effects of the small Am^/Amjj 
and sin #13 become important in the MSW resonance region. 

Finally, we describe two problems that we could not fully address in this paper, and which are tasks for 
future research. 

1. The approximate formulas in this paper are derived by using the condition #23 = 45°, which is the 
center value obtained from the atmospheric neutrino experiments. However, but differences from this 
value may exist within 90% confidence level. 

2. The condition for the symmetric matter density is satisfied in the 1-dimensional models, like the 
PREM and the akl35f, but the actual matter density, for example, that from J-PARC to Beijing is 
not symmetric |52j . Therefore, our aim for future work is, to derive more sophisticated approximate 
formulas that hold not only in symmetric matter but in arbitrary matter as well. 

To solve the above two problems are the future works This is now included in the upper sentence. 
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A General Feature of CP Dependence 

In this appendix we calculate the coefficients of the probabilities in detail. We show that the 2-3 mixing 
angle and the CP phase are not affected by matter, from a different point of view as described in our 
previous paper |41| . This result means that we only have to consider the matter effects on four parameters 
(Am|x, 612) and (Am^, 6*13)- By using this result, we can understand the matter effects in three generations, 
which become complex compared with that in two generations. 

A.l Remarkable Features of Effective Masses 

In this subsection, we show that (#23, $) do not affect the effective mass in three generation Hamiltonian. 
If we express the effective Hamiltonian in matter as 

H = C/diag(0, A21, A 3 i)^ + diag(a, 0, 0), (A.l) 

the equation of eigenvalue is given by 

det(t -H) = t 3 - (A 2 i + A 3 i + a)t 2 

+ (A21A31 + o(A a i(l - \U e2 \ 2 ) + A 3 i(l - \U e3 \ 2 )))t - aA 21 A 31 \U el \ 2 = 0, (A.2) 

and by solving this equation, we obtain the effective masses as 

Ai = ^ - \^ A2 - % BS - ^tV^ 2 - 35 y/l - S 2 (A.3) 
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A 2 
A, 



A 



3 3 v 



3BS 



|46| . where A, B, C and 5 are defined by 

A = A 2 i + A 3 i 



B = A21A31 + o[A a i(l - 1 17, 



<-2 



A3l(l-|C/e3| 2 )] 



C 



oAaiAailt/eil 



S = cos 



1 

- arccos 
3 



2A 3 - 9AB + 27C 



2yf{A 2 - 3B) 3 

These effective masses depend only on the following three vacuum mixing angles 

I Cell = C12C13, |Ue2| = S12C13, |?7 e 3|=Si3. 



(A.4) 
(A.5) 



(A.6) 
(A.7) 
(A.8) 

(A.9) 



(A.10) 



One can see from these equalities that the effective masses are independent of the 2-3 mixing angle #23 and 
the CP phase 5. Next, let us consider this result from a different point of view. 



A. 2 Decomposition of 2-3 mixing and CP Phase from Hamiltonian 

In this section, we separate #23 an d 6 from the Hamiltonian and we study the dependence of the 
amplitudes on the two small parameters a and S13. The Standard Parametrization is defined by 

c/ = 02 3 rOi 3 r t Oi2, (A.n) 

where the CP phase matrix T is given by 

r = diag(l, 1, e iS ). (A.12) 

The CP phase matrix T and the 1-2 mixing matrix On are commutable as 

[T,0 12 ] = [r,diag(0,A 21 ,A 3 i)] = 0. (A.13) 

Therefore, the Hamiltonian can be separated as 

H(t) = (7diag(0, A21, AsOtft + diag(a(f), 0, 0) = ChaTH'^O^, (A.14) 

where H'(t) is defined by 

H'(t) - Oi 3 Oi 2 diag(0, A21, A 3 i)Of 2 Of 3 + diag(o(t), 0, 0). (A.15) 

This means that the 2-3 mixing and the CP phase can be separated from the part which includes the matter 
effects a(t). 

In the case of constant matter density profile, we obtain 

det(A — H) = det(A — H'), (A.16) 

the 2-3 mixing angle and the CP phase do not affect the eigenvalue equation. Accordingly, the effective 
masses are independent of the 2-3 mixing angle and the CP phase, which coincide with the result obtained 
in the previous subsection. 



A. 3 Exact CP and 2-3 mixing Dependence of Oscillation Probabilities 

Here, let us consider the case in which we apply the above discussion used in the Hamiltonian to the 
amplitude. Solving the Schrodinger eq. for the amplitude in matter, we obtain 

S(t) = Tcxp j-i J H(t)dt^ . (A.17) 



22 



By using this, we obtain 



S(t) = Tcxp j-i J 23 TH'{t)T^0^ 3 dt\ = 23 rTexp j-i J H'(t)dt^ r+Of 3 = C^srS'ftjrtO^ (A.18) 

5(t) = 02srs'(t)r+o^. (A.i9) 



from QA.14J1 . Therefore, S(t) satisfies 
From this equation, we obtain 



P(v e ^v e ) = C ee , (A.20) 

P{v a -fv l} ) = A Q/3 cos<5 + B aP sm5 + C a/3 , (A.21) 

when the initial or final state is z/ e , and 

P(y a — > i/jg ) = cos (5 + _B Q /3 sin 5 + C Q) g + D a p cos 25 + £ aj g sin 2(5, (A. 22) 

in the case of v a ,vp — v^, v r |41| . The final result is given by 

P(v^v e ) = C ee = \S' e f, (A.23) 

P{ve^v^) = A efl co8 8 + B efl ain5 + C efl , (A.24) 

A e „ = 2Re[5^S; e ]pa3«23, (A.25) 

B eli = -2lm[S* e S' Te ]c 23 s 23 , (A.26) 

Ce» = \3^\ 2 4 3 + \S' T f4 3 , (A.27) 

P(v e ^v T ) = A eT cos 5 + B eT sin 5 + C eT , (A. 28) 

A eT = -2Re[^*5; e ]c 23 S23, (A.29) 

B eT = 2Im[^* e S; e ]c 2 3S23, (A.30) 

Cer - Kfsl 3 + \S' r fc 2 23 , (A.31) 

P(i^ — > = A w cos (5 + B MM sin (5 + C w + cos 25 + sin 2(5, (A.32) 

A m = 2M(S'^4s + S! rT s 2 23 r(S' Tll + S , flT )}c 2 ss 23 , (A.33) 

= -2Im[(S^cl 3 + S' TT sl 3 )^S^-S^)}c 23S23 , (A.34) 

C W = l^| 2 4 + l^r| 2 4 + (l^r| 2 + l^| 2 +2Rc[5;;^ r ])4 S 2 3 , (A.35) 

D w = 2Re[S;;s; T ]44 !J (A.36) 

= 2Im[5;;^ r ] C | 3S | 3 , (A.37) 

P(v T — ► f r ) = ^tt cos (5 + -B Tr sin (5 + C TT + D TT cos 2(5 + E TT sin 2(5, (A. 38) 

A TT = -2Rc[(^4 + ^ T ^ 3 )*(^ + ^ r )]c23 S 23, (A.39) 

b tt = 2im[(s; M 4 + s; T c| 3 )*(s; M -s;j]c23 S 23, (a.4o) 

Crr - |^| 2 4 + |5; r | 2 c4 3 + (|^ r | 2 + |^| 2 + 2Re[5;;5; r ])c2 3S 2 3 , (A.41) 

D TT = 2Rc[^;^ r ]c 2 3S 2 3 , (A.42) 

E TT = 2Im[5;;^ r ]c 2 3S 2 3 , (A.43) 

P(Vft — ► f r ) = ^4 M r cos (5 + S MT sin (5 + C MT + Z? Mr cos 25 + E^ T sin 25, (A.44) 

A, T = -2Rc[(S^-S' TT )%S' T ^ 3 -S^s 2 23 )}c 23 s 23: (A.45) 

B, lT = 21m[(^-^ r )*(^c 2 3 + ^ rS 2 3 )] C2 3 S 23, (A.46) 

C, r = |^ r | 2 4 3 + |^| 2 C 4 3 + (|^J 2 + |^ T | 2 -2Rc[5;;5; r ])c 2 3S 2 3 , (A.47) 
IV = -2Rc[^;^ T ]c 2 3S 2 3 , (A.48) 
E, T = -2lm[S'*S' }c 2 23 s 2 23 . (A.49) 
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